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ABSTRACT 



A variety of physical, social and biological systems generate complex fluctu- 
ations with correlations across multiple time scales. In physiologic systems, 
these long-range correlations are altered with disease and aging. Such corre- 
lated fluctuations in living systems have been attributed to the interaction of 
multiple control systems; however, the mechanisms underlying this behavior re- 
main unknown. Here, we show that a number of distinct classes of dynamical 
behaviors, including correlated fluctuations characterized by l//-scaling of their 
power spectra, can emerge in networks of simple signaling units. We find that 
under general conditions, complex dynamics can be generated by systems fulfill- 
ing two requirements: i) a "small-world" topology and ii) the presence of noise. 
Our findings support two notable conclusions: first, complex physiologic-like 
signals can be modeled with a minimal set of components; and second, systems 
fulfilling conditions (i) and (ii) are robust to some degree of degradation, i.e., 
they will still be able to generate l//-dynamics. 



3 



Complex systems are typically comprised of a number of interacting units that communicate in- 
formation and are able to process and withstand a broad range of stresses Jlll^lslQ]. In physiol- 
ogy, free-running healthy systems typically generate complex output signals that have long-range 
correlations — i.e., a 1 //-decay of the power spectra for low frequencies 1 J^QQ]. Deviations 
from the 1/f pattern have been associated with disease or aging on a number of contexts J^Hsl]. 

In spite of its practical and fundamental interest [9], the origin of such correlated dynamics 
remains an unsolved problem ||4j]. Until recently, attention has focused primarily on the complex- 
ity of the specific physiologic sub-systems or on the nature of the nonlinear interactions between 
them 111(1 IllL Il2ll . In particular, Boolean variables, which can take one of two values, or 1, 
and Boolean functions have been extensively used to model the state and dynamics of complex 
systems — see I12I for an introduction. The reason such a "simplistic" description may be ap- 
propriate arises from the fact that Boolean variables provide good approximations to the nonlin- 
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ear functions encountered in many control systems BIOL I13LI14L 11511 . Random Boolean networks 
(RBNs) were proposed by Kauffman [KJ as models of genetic regulatory networks, and have also 



been studied in a number of other contexts [13, LLJ]. Wolfram [15], in contrast, proposed that 
cellular automata (CA) models — a class of ordered Boolean networks with identical units — may 
explain the real-world's complexity. Neither of these two classes of models has been shown to 
generate the complex dynamics with l//-fluctuations observed in healthy physiologic systems. 
We propose here a new modeling approach (Fig. la) that departs from traditional approaches 



in that we pay special attention to the topology of the network of interactions [4], and the role of 
noise Jl^l . Our model is rooted in two considerations frequently observed in real- world systems: 



(i) the units in the system are connected mostly 



ocally but also with some long-range connections, 



giving rise to so-called small-world topology H 1 7L 1 1 811 . and (ii) the interaction between the units 



2311 . We demonstrate 



is affected by noisy communication and/or by noisy stimuli [19, 
that simple rules, such as the majority rule, are able to generate signal with complex fluctuations 
under simple but physiologically-relevant conditions. 

Methods 



1 A power law decaying power spectrum S(f) oc f^ 13 is the signature of a signal with power-law decaying 
correlations. The case j3 = 2 corresponds to a Brownian, while (3 — corresponds to a completely 
uncorrelated "white" noise. The intermediate case, S(f) oc l/f, is a "compromise" between the small 
time-scale roughness but large time-scale smoothness of white noise and the small time-scale smoothness 
but large time-scale roughness of Brownian noise. 
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FIG. 1: Emergence of complex dynamics in simple signaling networks. Our model for signaling networks is 
defined as follows, a, The units comprising the network, which are located on the nodes of a one-dimensional 
lattice, have bi-directional nearest-neighbour connections, b, A number k e N of additional unidirectional 
links is established between pairs of randomly selected units, where k e is the mean excess connectivity and 
N is the number of units in the system. At time t = 0, we assign to each unit i — 1, . . . , N a state <Tj(0) 
randomly chosen from the set {0, 1} and a Boolean function Ti\ cf. Eq. (1). This Boolean function (or rule) 
determines the way the inputs are processed. Each unit effectively processes two inputs, one corresponding 
to the average state of its neighbours and its own state. With probability 77, a unit "reads" a random 
Boolean variable instead of the state of a neighbour, where the parameter r\ quantifies the intensity of the 
noise. Note that the noise does not alter the state of the units, just the value read by its neighbour. At each 
subsequent time step, each unit updates its state synchronously according to its Boolean function. Time 
evolution of systems comprising 512 units with T% = 232 for all units and 77 = 0.1, and c k e — 0.15, d 
k e — 0.45, e k e = 0.90. A red dot indicates o^(i) = 1 while a yellow dot indicates <Zi(i) = 0. The three 
panels display the time evolution for systems starting from the same initial configuration and using the same 
sequence of random numbers. Thus, the difference in the dynamics is due uniquely to the different number 
of long distance links. For k e = 0.15 the system quickly evolves toward a configuration with several clusters 
in which all the units are in the same state. The boundaries of these clusters drift due to the noise but 
the state of the system S(t) is quite stable and the dynamics are close to Brownian noise. In contrast, for 
k e — 0.90 a large stable cluster develops and the state of the system changes only when some units change 
state due to the effect of the noise. This process yields white-noise dynamics. For k e = 0.45 clusters are 
formed but they are no longer stable — unlike what happens for small k e . In this case, information propagates 
through the random links which can lead to a change in the state of one or more units inside a cluster. Our 
results suggest that because these long-range connections exist on all length scales, they lead to long-range 
correlations in the dynamics and the observed l//-behavior (Fig. 3b). 

The model — We place the Boolean units comprising the network on the nodes of a one- 
dimensional ring and establish bi-directional nearest-neighbor connections (Fig. 1). Then, we 
add k e N additional unidirectional links — where k e is the mean excess connectivity and N is the 
number of units in the system — between pairs of randomly selected units. Hence, each unit has 
a set of links through which incoming signals arrive and that the unit then processes. The state 
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256 Boolean rules 
of 3 inputs 



64 symmetric rules 

generalizable to 
any number of inputs 
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FIG. 2: Selection of Boolean rules for investigation. Our goal is to investigate Boolean functions that 
display non-trivial dynamics and can be generalized to any number of inputs. To this end, we start from 
256 rules of 3 inputs but then restrict our attention to the ones that are symmetric under permutations of 
the external inputs. This selection results in 64 Boolean rules. However, each rule has another that is its 
complementary — i.e., display the same dynamics when switching zeros and ones — or inverse — i.e., display 
the same dynamics when taken every other step. Since these pairs of rules have equivalent dynamics, we need 
to investigate only 32 independent rules. Of these, eight do not display fluctuations, even in the presence of 
noise, resulting in 24 independent rules that could present complex fluctuations. The phase-spaces of each 
of these 24 rules are shown in the Supplementary Material. 

of each neighbour is replaced by a random value with probability r\ — which parameterizes the 
intensity of the noise. 

We assign to each unit i = 1, . . . , N a Boolean function jF i5 which determines the way the states 
of the neighbors and its own state are processed. We restrict our study to Boolean functions that 
have only two "effective" inputs: the state of the unit and the average state of all other neighbours. 
This restriction yields 64 unique symmetric Boolean functions (see Fig. 2 and the Supplementary 
Material) and it has the advantage that it permits a topology-independent implementation of the 
Boolean functions, thus enabling a systematic study of the effect of different rules on the dynamics 
of the system. 

Quantification of the dynamical behavior of the system — We start all of our numerical simu- 
lations with a random initial configuration and let the system evolve synchronously according to 
the rules of the model. We define the state S(t) of the system as the sum of the states o~i of all the 
Boolean units 

S(f)=5>(f). (1) 

% 

We record S(i) during the course of the simulation (see Figs. 3a-c) and quantify the complexity 
of the time series generated in terms of its auto-correlation function [3]. We apply the detrended 



LQJ UJ L2J L3J 



Select ruleslhat are" 
symmetric under 
permutations of 
external inputs 



Jlemove 
complementary rules 




CE CJJ HD (X) 



6 



3 

d 

tu 

00 



I 1 I 1 < 1 < 1 
. , . i . i . e i 




i 1 i 






1 , | . | , | , 

h-l™ & =0.45 


i | i | i | i 

b2 


1 1 
i 


c 


1 ' 1 
1 1 


1 1 1 
. i . e 


1 1 

= 0.90 




1000 2000 3000 4000 5000 

Time step 



10 10~ 

Time scale, n 



FIG. 3: Quantification of the correlations in the state of Boolean signaling networks. As discussed in the 
text, we define the state of the system as S(t) — Yl a i(t)- We show S(t) for a system with N = 4096 
units, rj = 0.1, Ti = 232, and k e — 0.90, k e = 0.45, k e = 0.15. The three values of k e lead to quite 
different dynamics of the system, a, For a small number of random links, the time correlations display 
trivial long-range correlations such as found for Brownian noise, b, For an intermediate value of k e , long- 
range correlations emerge and the power spectrum displays a power-law behavior, S(f) oc l//* 3 , with /3 ml. 
Panels bl-b2 display the state of the system according to different definitions: In bl the state of the system 
is defined as the sum of the states of a random sample comprising l/8th of all units, whereas in b2 the 
state of the system is defined as the sum of the states of a block of contiguous units comprising 1 / 8th of the 
systems. Our results indicate that the evolution of a subset of the population is similar to the dynamics of 
the whole system, c, For a large number of random links, k e = 0.90, the dynamics are less correlated, d, 
Estimation of temporal auto-correlations of the state of the system by the detrended fluctuation analysis 
(DFA) method [5j. We show the log-log plot of the fluctuations F(n) in the state of the system, versus time 
scale n for the time series shown in Figs. 3a-c. In such a plot a straight line indicates a power-law dependence 
F(n) oc n a . The slope of the lines yields the scaling exponent a, which for a number of physiologic signals 
from free-running healthy mature systems take values close to one [3j. The exponent a is related to the 
exponent f3 of the power spectrum of the fluctuations, S(f) oc l//* 3 , through the relation (3 = 2a — 1. The 
data sets have been shifted upward and the different sets correspond, from top to bottom, to the time series 
shown in Panels 2a-c. 

fluctuation analysis (DFA) method [5] which quantifies long-range time-correlations in the dy- 
namical output of a system by means of a single scaling exponent a (Fig. 3d): Brownian noise 
yields a = 1.5, while uncorrelated white-noise yields a = 0.5. For many physiologic signals, one 
observes a ~ 1, corresponding to 1 //-behavior, which can be seen as a "trade-off between the 
two previous cases fl. 

Results 

Random Boolean networks — The RBN model corresponds to a completely random network with 
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randomly selected Boolean rules for the units. As shown in Fig. 4a, we find white-noise dynamics 
for essentially any pair of values of k e and r\ within the ranges considered, suggesting that even 
in the presence of noise a system of random Boolean functions cannot generate l//-dynamics. 
This result is not unexpected, since the random collection of Boolean functions comprising the 
system prevents the development of any order or predictability in the dynamics. 

CA models with small-world topology and noise — We systematically study the 64 symmetric 
rules (see Supplementary Material and Fig. 2) for different pairs of values of k e and 77. Some of 
the rules have parallels to physiologically meaningful dynamics. Rule 232 is a majority rule, that 
is, each unit will be active next time-step only if the majority of its neighbors is active now. Rule 50 
is a threshold rule with refractory time period, that is, whenever the inputs of the neighbors surpass 
a certain value a unit becames active in the next time-step and then will be inactive for at least one 
time-step. The 64 symmetric rules lead to three qualitatively distinct phase-spaces (Figs. 4b-d). 
Rule 232, the majority rule — which is representative of the first type of phase-space — displays 
three distinct types of dynamical behavior (Fig. 4b): For small k e , we find mostly Brownian-like 
scaling. For large k e , we find mostly white-noise dynamics. Of greatest interest, for intermediate 
values of k e and for a broad range of values of the intensity of the noise 77, we find 1 / /-fluctuations. 

Rule 50 — which is representative of a second type of phase-space — displays fewer types of 
dynamical behavior. In particular, we find only a narrow range of noise intensities (with a weak 
dependence on k e ) for which the dynamics display 1/ /-correlations. For n ~ 0.1, the dynamics 
become uncorrected (Fig. 4c). Rule 160 — which is representative the third type of phase- space — 
displays white-noise dynamics for all values of k e and 77 (Fig 4d). 

Note that for k e = — i.e., when the network is a one-dimensional lattice — the model is not able 
to generate l//-dynamics. This implies that in the context of the model, the existence of long- 
range connections — i.e., the small-world topology achieved by making k e 7^ — is an essential 
ingredient for the emergence of 1/ /-dynamics. 

Robustness of the findings 

In order to determine the generality of the results presented above, one needs to address the ques- 
tions of how these findings are affected by (i) changes in the topology of the network or (ii) "errors" 
in the units' implementation of the rules. 

Concerning (i), we note that the network topologies considered so far span the cases of ordered 
one-dimensional lattices, small- world networks, and random graphs [24]. However, all networks 
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FIG. 4: Systematic evaluation of the correlations in the dynamics generated by different rules. We quantify 
the long-range correlations in the dynamics by means of the DFA exponent a |5j systematically estimated for 
time-scales 40 < n < 4000. We show a for 3,721 pairs of values of k e and the noise r\ in the communication 
between the units comprising the network. For all simulations, we follow the time evolution of systems 
comprising 4,096 units for a transient period lasting 8,192 time steps, and then record the time evolution 
of the system for an additional 10,000 time steps. In order to avoid artifacts due to the fact that for some 
of the rules the units switch states with period 2, we consider in our analysis the state of the systems at 
every other time step, a, Random Boolean network (RBN) as defined by Kauffman 10]. Our results show 
that the dynamics generated by these systems are generally of the white-noise type with a weak dependence 
on the noise intensity and no dependence on the number of long-distance links, b, Rule 232, a.k.a. the 
majority rule. This rule is representative of two other rules, numbers 19 and 1. Rule 232 displays a very 
rich phase-space with a variety of dynamical behaviors, all the way from white noise (represented in white 
and green) to Brownian noise (represented in black), c, Rule 50, a threshold rule with refractory period. 
This rule is representative of eight other rules, numbers 5, 36, 37, 73, 77, 94, 108, 164. These rules display 
a relatively simple phase-space with behaviors extending from white noise to l//-noise. The l//-behavior 
is restricted to very small noise intensities and there is a very weak dependence on k e . c, Rule 104. This 
rule is representative of twelve other rules(see the Supplementary Material). Their phase space is extremely 
simple, displaying only white noise behavior. 

considered are comprised of units with approximately the same degree, i.e., the same number of 
connections. To investigate the role of the distribution of number of connections, we also study 
networks which span the range of empirically observed degree distribution: a delta-distribution, an 
exponential distribution, and a power law distribution. The latter case corresponds to the so-called 
scale-free networks HX8I1 - Notably, we find that the picture of the phase-space presented in Fig. 4b 
does not get altered by these changes in the degree distribution (see the Supplementary Material). 

In order to address (ii), we systematically explore the dynamical behaviors in the phase-space 
defined by (k e , 77) for systems composed of units operating according to either rule 232 or a ran- 
domly selected rule. Notably, we find that with as many as l/4th of all units operating according 
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FIG. 5: Phase-space for signaling networks with mixing of Boolean rules. We systematically calculate 
the exponent a, characterizing the correlations in the dynamics, for systems composed of units operating 
according to rule 232 but with some fixed fraction of units operating according to a randomly selected 
symmetric Boolean rule. Each value is an average over 5 independent runs, a, l/16th of the units operating 
according to a randomly selected rule, b, l/8th of the units operating according to a randomly selected 
rule, c, l/4th of the units operating according to a randomly selected rule, d, l/2th of the units operating 
according to a randomly selected rule. These figures suggest that the presence of random Boolean functions 
leads to a decrease in the richness of the phase space of the systems. Specifically, if more than l/4th of all 
the units operate according to randomly selected Boolean functions, then the phase-space displays mostly 
white-noise dynamics. 

to random Boolean functions the model still displays a rich phase-space including white, 1/ / and 
Brownian noise (Fig. 5). This finding holds even if instead of using random Boolean functions we 
consider a single Boolean function (Fig. 6). 

Discussion 

Our results are notable for a number of reasons. First, they demonstrate that a new model of 
signaling networks is able to generate a broad range of behaviors reminiscent of those observed 
in physiologic systems. Second, we show for a rather general class of models that a variety of 
dynamical behaviors can only emerge under restrictive, but physiologically-relevant assumptions; 
namely, the system must have a small- world topology, and noise must be present. 

An interesting aspect of our results is that some of the rules we consider have plausible phys- 
iologic interpretation. For instance, rule 232 is a rule in which a unit changes its state to that of 
the majority of the incoming inputs. A majority rule appears to be operative in the central nervous 
system where multiple fibres (excitatory or inhibitory) converge onto a single neuron. Action po- 
tentials converging in a neuron summate to bring the neuron to the threshold for firing [25D. A 
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FIG. 6: Phase-space for signalling networks with mixing of two Boolean rules. We systematically calculate 
the exponent a, characterizing the correlations in the dynamics, for systems composed of units operating 
according to either rule 232 or rule 50. Each value is an average over 5 independent runs, a, l/16th of 
the units operating according to rule 50 and 15/16th operating according to rule 232. c, l/8th of the units 
operating according to rule 50 and 7/8th operating according to rule 232. d, l/4th of the units operating 
according to rule 50 and 3/4th operating according to rule 232. e, 1/2 of the units operating according to 
rule 50 and 1/2 operating according to rule 232. When both rules are present in the system — and at least 
50% of the units operate according to rule 232 — we still find several distinct classes of dynamical behaviours, 
including a wide range of parameter values that generate 1//- noise. 



majority rule also appears to be operative in the baroreflex control of the cardiovascular system. 
The baroreflex is a feedback loop that continuously controls heart rate by modulating the degree 
of sympathetic and parasympathetic nervous system input to the sinus node of the heart. Changes 
in heart rate values are determined by whichever input is dominant at the moment [25]. 

Our results demonstrates that complex fluctuations are present even when a fraction of the units 
obeys randomly chosen Boolean functions. This suggest that systems comprised of units operating 
according to the majority will be robust to the removal or failure of units. Additionally, since 
this complex dynamics does not depend on a scale-free topology, this system do not necessarily 
displays the vulnerability to targeted attack observed in scale-free networks J3]. In biological 
systems this robustness could support the "physiologic reserve" enabling an organism to overcome 
age- or disease-related loss of system components. 

Additionally, our findings raise the intriguing possibility that the interactions within physiologic 
systems and their degradation with aging and pathology may be symbolically mapped as a "walk" 
on the phase space (k e ,r]). According to this model, 1 / /-dynamics similar to that found for healthy 
physiology are generated when noise intensity and connectivity reside in a well-defined range J^ . 
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A loss of complexity with a breakdown of long-range correlations would be anticipated when these 
parameters assume values outside this range (Fig. 4). Support for this formulation comes from 
analysis of heart rate dynamics with aging and disease, where connectivity or coupling among 
system components is likely to be degraded J^l 3, Similarly, evidence suggests that decreased 
social connectedness, and the corresponding decrease in "noisy" stimulation, may be associated 



with increased cardiac mortality and decreased functional recovery from stroke or dementia [ 28 

Q. 




In a related way, the development of complex, adaptive dynamics during the maturation of the 
organism may be accounted for, at least in part, by the evolution of appropriate connections (see 
I30II31I for empirical evidence). Our model predicts that the fraction of non-local connections has 
an optimal range of values; hence, an excessive number of certain types of inputs may also degrade 
functionality. Our results are thus consistent with empirical evidence suggesting that development 
and maintenance of healthy function may require adjusting the number of connections. Finally, the 
model may also provide a robust way to generate fluctuations that closely resemble physio 
signals, which could be implemented in medical devices such as mechanical ventilators 
and hormone infusion pumps II34I1 . 
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